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Abstract 

The effects of a temperature-dependent rheology on large-scale continental 
extension are investigated using a thin viscous sheet model. A vertically- averaged 
rheology is used that is consistent ■with laboratory experiments on power-law 
creep of olivine and that depends exponentially on temperature. Results of the 
calculations depend principally on two parameters: the Peclet number, which 
describes the relative rates of advection and diffusion of heat, and a dimensionless 
activation energy, which controls the temperature dependence of the rheology. 
At short times following the beginning of extension, deformation occurs with 
negligible change in temperature, so that only small changes in lithospheric 
strength occur due to attenuation of the lithosphere. However, after a certain 
critical time interval, thermal diffusion lowers temperatures in the lithosphere, 
strongly increasing lithospheric strength and slowing the rate of extension. This 
critical time depends principally on the Peclet number and is short compared 
with the thermal time constant of the lithosphere. The strength changes cause 
the locus of high extensional strain rates to shift with time from regions of high 
strain to regions of low strain. Results of the calculations are compared with 
observations from the Aegean, where maximum extensional strains are found in 
the south, near Crete, but maximum present-day strain rates are largest about 
300 km further north. The comparisons support the hypothesis that the observed 
separation of the loci of maximum strain and maximum present-day strain rates 
in the Aegean may be the consequence of changes in lithospheric strength arising 
from the temperature-dependent mechanical properties of lithospheric materials. 
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1. Introduction 

Recently, much attention has been directed towards understanding con- 
tinental extensional tectonics. Simple kinematic models of stretching continental 
lithosphere (e.g., McKenzie, 1978a; Jarvis and McKenzie, 1980; Royden and Keen, 
1980) have been shown to be consistent with the subsidence histories of a number 
of sedimentary basins and passive margins (e.g. Sclater and Christie, 1980; Roy- 
den and Keen, 1980; Le Pichon and Sibuet, 1981; Royden et al., 1983). In these 
models, horizontal stretching of the lithosphere (and attendant vertical thinning) 
perturbs the geotherm by raising the thermal gradient and by introducing hot 
asthenospheric material from below to replace attenuated lithosphere. The ther- 
mal evolution of the lithosphere as it stretches depends principally on the relative 
rates of advection of heat into the lithosphere by upwelling of asthenosphere and 
diffusive loss of heat through increased surface heat flow, while thermal diffusion 
dominates after the end of stretching. 

A logical extension of thermal models of lithospheric extension is the 
consideration of dynamic effects of rheology on the deformation. Earth materials 
generally have temperature-dependent material properties, so the thermal evolu- 
tion and the mechanical evolution of stretching lithosphere are strongly coupled. 
In general, diffusion of heat during stretching results in the cooling of individual 
horizons in the lithosphere, as they are brought close to the earth’s surface; with 
continued stretching this cooling may cause the vertically-averaged strength of 
the lithosphere to increase (England, 1983). This suggests that there may be a 
limit to the extensional strain attainable which would depend primarily on the 
magnitude of the driving stresses, the rate of deformation, and the degree of tem- 
perature dependence of the rheology (Houseman and England, 1986a; Sonder et 
al., 1987). 
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The calculations summarized above all considered the evolution of a sin- 
gle vertical column of lithosphere, and for simplicity, consideration was not given 
to the thermal and mechanical evolution, in three dimensions, of extending con- 
tinental lithosphere or the effects thereon of a temperature-dependent rheology. 
This paper represents a preliminary investigation of the effects of such a rheology 
on the three-dimensional deformation of the continents in an extensional enviro- 
ment. We use a thin viscous sheet model (Bird and Piper, 1980; Vilotte et al., 
1982, 1986; England and McKenzie, 1982, 1983) and treat the lithosphere as a 
continuum with vertically averaged properties. We first present general results of 
the calculations that illustrate the coupling between the thermal and mechanical 
behavior of the thin sheet, and we identify the phenomena that arise due to the 
inclusion of the temperature dependence. 

The most important results presented in this paper concern the distribution 
of strains and strain rates. Calculations that include a temperature-dependent 
rheology result in a separation of the loci of maximum strain and maximum 
strain rate. The separation occurs as a result of increasing mechanical strength in 
thinned areas and begins after a time that can be approximated by the time 
required for thermal conduction to first affect temperatures in the uppermost 
mantle. Calculations with little or no temperature dependence do not produce 
this pattern of deformation. The results are compared with observations from 
the Aegean that indicate that, while the maximum Neogene extensional strain 
occurred in the southern Aegean north of Crete (Angelier et al., 1982), the max- 
imum present-day strain rates are found in the northern Aegean. The consistency 
of the results of calculations with the observations allows us to propose a simple 
mechanism by which this deformation pattern may have been produced. 
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2. Choice of rheology 

The thin viscous sheet model used in this paper requires the use of a verti- 
cally averaged rheology, and the choice of such a rheology, which is also to 
include temperature dependence, is of primary importance to the results dis- 
cussed below. Laboratory experiments on olivine (e.g. Goetze, 1978) suggest that 
much of the long-term strength of the continental lithosphere may reside in the 
uppermost mantle (see Brace and Kohlstedtj 1980). Therefore, we are motivated 
to choose a temperature-dependent rheology based on olivine, and will neglect 
the strength of the crust. Although this is necessarily an approximation to the 
behavior of geological materials over the range of temperatures, pressures, strain 
rates, and compositions relevant to continental deformation, the adoption of a 
more realistic, but more complicated, rheology is unwarranted at this stage. Use 
of the simple rheology based on the physical properties of olivine avoids the need 
to introduce extra parameters describing the dependence of the strength of the 
crust on temperature and, especially, the nature of the poorly understood mid- 
crustal transition from frictional sliding to creeping flow. 

The laboratory results for olivine suggest an empirical creep law of the form 

k = A expj-^j (g x -< 7 3 ) n (1) 

(e.g. Goetze, 1978) for differential stress g x -cr 3 less than 200 MPa. Q is the 
activation energy, R is the gas constant, T is absolute temperature, e is strain 
rate, and A is a constant. The stress exponent n typically has a value of about 

3. Generalizing in terms of tensor components of strain rate and deviatoric stress 
(r), and writing as a vertical average, Eq. 1 becomes 



where E is the second invariant of the strain rate tensor, equal to (e,y e,y ) 1//2 
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(with summation over repeated subscripts), and r,y and c,-y represent vertically 
averaged quantities. B includes the temperature dependence and will be 
referred to below as the strength of the lithosphere: 

B = c / exp ijrfc) iz (3) 

where z L and z M are the depths to the base of the lithosphere and the Moho, 
respectively, and C is a constant. 

In a later section we present results of calculations of Moho temperature, 
The significance of the Moho temperature is that it is a convenient indicator 
of lithospheric strength for continental lithosphere whose deformation is governed 
by a power law such as Eq. 2 (England, 1983; Sonder and England, 1986). This 
may be seen by considering an approximation to Eq. 3. For a constant thermal 
gradient 7, the integral in Eq. 3 may be evaluated: 



England (1983). Thus, the strength increases sharply with decreasing Moho tem- 
perature. To a much lesser extent, B depends on the thermal gradient, <7; this is 
important if thermal diffusion is small compared with advection of heat, when 
deformation may alter the thermal gradient without appreciably changing the 
Moho temperature. 

A second important parameter is Q /RT^, where T L is a reference tempera- 
ture, here taken as the temperature at the base of the lithosphere. For the magni- 
tude of Q suggested by laboratory experiments on olivine, Q /RT L falls in the 
range 30-50, if we assume that T L is around 1400 0 C. 
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3. Duration of extension: simple analytic model 

England (1983) showed that, during continental extension, decreasing tem- 
peratures in strength-controlling regions of the lithosphere (e.g. the uppermost 
mantle) can cause the vertically-averaged strength of the lithosphere to increase 
sharply. For a range of conditions of driving stress, strain rate, and thermal state, 
this strength increase can be sufficient to limit further extension (e.g., Houseman 
and England, 1986a). 

If extension of the lithosphere were to occur instantaneously and homogene- 
ously by a factor /?, there would be a discontinuity in the slope of the geotherm 
at the base of the lithosphere (depth z L /0), which would relax with time, result- 
ing in cooling of the lithosphere (see Figure 1) (McKenzie, 1978a). The exponen- 
tial dependence of strength B in Eq. 4 on the temperature at the Moho means 
that rather small changes in Moho temperature produce large changes in B (Eng- 
land, 1983; Houseman and England, 1986a). Therefore, to a good approximation 
we can neglect changes in temperature except as they appear in this exponential 
term and write that the change in Moho temperature AT (6 ) required to produce 
a change in lithospheric strength by a factor b is 

tiR T 

A 7-( t )«__iln(6) (5) 

Evaluation of Eq. 5 for values of the parameters believed to be appropriate to the 
continental lithosphere (Qw 5X10 2 kJ mol -1 ; n = 3; T M « 800 K shows that a 
ten-fold increase in lithospheric strength (6 = 10) results from a drop of a few 
percent in Moho temperature, so the approximations in going from Eq. 4 to Eq. 5 
are not severe. 

Because we are concerned with small temperature changes at the Moho, and 
hence with times small compared with the time for thermal relaxation of the 
lithosphere, the time for these temperature changes is governed by the distance 
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( z l~ 2 m)/P from Moho to the "knee” in the thermal profile (Figure 1). 
Neglecting geometrical factors, the change in Moho temperature may be written: 


AT 


Ti erfc 


( Z L ~ Z M ) /£ 

2n//cF 


( 6 ) 


where k is the thermal diffusivity. Eliminating A T between Eqs. 5 and 6 gives 
an expression for the length of time required to increase the strength of the litho- 
sphere by a factor b : 


( Z L ~ Z M ) /£ 

2 VkT 


erfc 1 



( 7 ) 


Evaluating Eq. 7 with Q ^ 5X10 5 J mol -1 K _1 , n = 3, and T i = 1600 K , and 
a Moho temperature approximately half that of the base of the lithosphere gives 

o 

the result that an increase in lithospheric strength by a factor of ten (6 = 10) is 
obtained at a time t c , where 

{*l~ z m) 2 /P 2 


Although this discussion relates to an instantaneously stretched lithosphere, 
it has the advantages of showing that the timescale for the evolution of strength 
of the lithosphere is generally short compared with the thermal time constant of 
the lithosphere as a whole, and of relating the timescale to the rheological param- 
eters. If, instead, we consider extension that occurs over a finite time interval, we 
may suggest a modification to the previous discussion. For simplicity, assume 
that the strain rate e is a constant equal to e 0 ; then 

0 = exp ( e 0 t) (9) 


and 
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<« (10) 

The form of Eq. 10 is a natural one for considering thermal diffusion alone, but a 
second timescale is involved when the mechanical problem is considered. In the 
calculations below, the extensional strain rates — at least before the development 
of rheological contrasts of the kind just discussed — depend on the characteristic 
horizontal length of the extensional boundary and on its speed (see England et 
al., 1985). For example, in the configuration we adopt here, strain rates near the 
extensional boundary are about v 0 /w, where t; 0 is the velocity of the boundary 
and tv is the width of the moving boundary (see Figure 3). In the mechanical cal- 
culations we use a dimensionless time t ' given by 


V 


v 0 t 

~D~ 


where D is the width of the solution domain. Eq. 10 may be written in dimen- 
sionless form as 


t 


t 

C 


Pe V 2 
8 




( 11 ) 


where lengths are made non-dimensional by = z /£> , and strain rates by 
e' — kD /v 0 . The Peclet number, Pe, in Eq. 11 is given by 



K 


( 12 ) 


It is this parameter that determines the relative importance of thermal diffusion 
and advection in the calculations below. At low Peclet number, thermal diffusion 
will influence the temperature structure, and hence the strength, before appreci- 
able strain has occured, whereas at high Peclet number, large strains may occur 
essentially isothermally. 
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Figure 2 shows solutions of Eq. 10 for values of 2, \/2, and 1 for the right 
hand side of Eq. 7, corresponding to values of 16, 8, and 4 for the constant in the 
denominator of Eq. 10. The right hand side of Eq. 7 can be easily evaluated for a 
given choice of values for the parameters. The three values of Q /RT L used in 
this paper are 50, 25, and 12.5 (Table 1). These, with b = 10 and the values of 
Tjv/ and Ti shown in Table 2 give values, respectively, of 1.5, 1.3, and 1.1 for 
the right hand side of Eq. 7; solutions of Eq. 10 using these values are illustrated 
in Figure 8, and compared with the equivalent times determined from three- 
dimensional calculations in which the changes in rheology explicitly affect the 
strain rate field. 

4. Three-dimensional calculations of extension with a temperature- 
dependent rheology 

4.1. Description of model 

Following England and McKenzie (1982, 1983), we approximate the deforma- 
tion of continental lithosphere by the deformation of a thin sheet of incompressi- 
ble viscous fluid that overlies an inviscid half-space. Because the horizontal 
dimensions of the sheet are much greater than the vertical dimension, and shear 
stresses are assumed to be zero on the top and bottom of the sheet, the deforma- 
tion can be approximated by neglecting vertical gradients of horizontal velocity 
components. This allows the deformation of the sheet to be expressed in terms of 
vertical averages of strain rate, deviatoric stress, and rheology. The sheet is also 
assumed to be in isostatic equilibrium. 

The problem consists of two parts, mechanical and thermal, that are coupled 
through the temperature dependence of the rheology. The formulation of the 
equations of motion, which follows England and McKenzie (1982, 1983), and of 



-.11 - 


the equation governing heat transfer are summarized in the Appendix. Solution 
of the equations is accomplished using a Lagrangean finite element scheme 
(Houseman and England, 1986b), also summarized briefly in the Appendix. 

4.2. Boundary and initial conditions and parameter ranges 

Figure 3 shows the boundary conditions used in the calculations. The solu- 
tion domain is initially a 2X1 rectangular box. Three sides {x' = -1, x' = 1, 
and y' = 1) are rigid ( u' = 0, v' =0). On the boundary initially at y' =0, 
the tangential velocity is zero; the normal velocity is specified such that part of 
the boundary moves outward with respect to the rest of the box, resulting in hor- 
izontal stretching and vertical thinning of the fluid near the boundary. 

Thermal boundary conditions are T' = 0 and T' =1 on the upper and 
lower surfaces, respectively, of the fluid. Initial conditions are uniform thickness 
of crust and lithosphere and a constant vertical temperature gradient at all 
points (x ,y ) in the box. Hence, initially, there are no horizontal gradients of layer 
thickness or temperature. 

Calculations were carried out using the values of the dimensionalizing con- 
stants in Table 1. For a given set of boundary conditions, there are four dimen- 
sionless parameters that define the solution space of the thin viscous sheet model 
with temperature-dependent rheology. Three of them have been introduced ear- 
lier (Q /RT l , Pe , and n ), and describe, respectively, the degree of temperature 
dependence of the rheology, the relative rates of advection and diffusion of heat, 
and the nonlinearity of the stress dependence of the rheology. Apparent values 
of these parameters appropriate to the deformation of the continents are dis- 
cussed more fully in the Appendix, and are listed in Table 2. We only consider 
deformation with n = 3 because this is the value appropriate to the power law 
creep regime -- in which the uppermost mantle is believed to be at geologically 
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relevant strain rates. The remaining dimensionless parameter is the Argand 
number ( Ar ; England and McKenzie, 1982; also see Section 4.3.4), "which 
expresses the degree to which buoyancy forces due to horizontal density gradients 
influence the deformation. 

4.3. Results 

We first present results with Ar = 0, in order to understand those aspects 
of the deformation that depend on the temperature dependence of the rheology. 
Calculations are presented later (Section 4.3.4) that have non-zero Argand 
numbers, so that the modifying effects of buoyancy forces on the deformation can 
be identified. 

4.3.1. Time evolution of selected nodes 

Figure 4 shows the time evolution of Moho temperature ( T M ), strength ( B ), 
vertical strain rate (e zz ), and strain (/?) for five nodes of the finite element mesh. 
For this experiment, Pe = 300, Q / RTi = 50, n == 3, and Ar = 0. Because of 
the Lagrangean finite element formulation used to solve the equations of motion, 
the mesh deforms with the fluid, so that by tracking these five nodes through 
time, we are tracking the history of five individual fluid parcels. 

At first, Moho temperatures and strengths are uniform over the entire solu- 
tion domain, so the distribution of strain rates is governed entirely by the boun- 
dary conditions and the value of the power law exponent (see England et al., 
1985). The region of greatest extensional strain rate is located on the boundary 
of the mesh, and the rate of extension decays with distance from the boundary. 
Extension proceeds for approximately 0.1 dimensionless time unit without a large 
change in Moho temperature (Fig. 4a). There is a slight (20%) decrease in 
strength (Fig. 4b), most apparent for the nodes closest to the extensional 
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boundary; this decrease is attributed to the increased thermal gradients resulting 
from the lithospheric thinning (see Eq. 4). Weaker areas, those closest to the 
boundary, concentrate the strain, so strain rates increase in those regions and 
decrease slightly elsewhere (Fig. 4c). 

With continuing extension, by t ' = 0.1 sufficient time has elapsed that the 
temperature at the Moho begins to drop (Fig. 4a). The temperature decrease is 
largest for nodes closer to the boundary, where there has been more extension. At 
the same time the strength begins to increase. For nodes 1, 2, and 3, the strength 
increase results in a sudden and sharp decrease in the rate of extension. Extension 
begins to be taken up in regions farther from the boundary (e.g. nodes 4 and 5) 
where the strength is less. By t ' = 0.25, the region of maximum rate of exten- 
sion is well within the interior of the fluid and the rate of extension on the boun- 
dary is approximately 10% of its original rate. 

As time proceeds, the strength also increases for nodes located away from 
the boundary (nodes 4 and 5), so strain rates peak there at t c* 0.2. The distri- 
bution of strain rates becomes spatially more uniform as the high strength of 
many nodes causes quasi-rigid behavior of the most highly extended parts of the 
fluid. By the end of the experiment, the strain rates for all 5 nodes are close to 
zero (Fig. 4c) and little additional strain is accumulating in the region illustrated 
in Fig. 4 close to the boundary (see Fig. 4d, in which /? has levelled off at values 
between 1.4 and 1.6). Almost all the deformation is being taken up further into 
the interior of the fluid. 

4.3.2. Strain rates 

In Figure 5 the horizontal distribution of vertical strain rate as a function of 
time is shown for the experiment of Figure 4 ( Pe = 300, Q /RT L = 50, 
n = 3, and Ar = 0). At first, extensional strain rates are greatest near the 
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boundary of the fluid; the form of the boundary condition also concentrates 
strain rates near the taper in the boundary velocity (Fig. 3). The early concen- 
tration of high strain rates near the boundary due to the slight strength decrease 
is plainly seen as a decrease in the spacing of contours. By t' = 0.2 the pattern 
of strain rate distribution is beginning to change. The rate of extension begins to 
drop, initially on the boundary where the strains are greatest, and the region 
with | e Z2 | > 0.8 begins to broaden and move towards the interior of the fluid. 
It is well away from the boundary by t' — 0.4, the end of the calculation. 

Figure 6 illustrates the influence of the choice of Peclet number and Q /RT L 
on the distribution with time of the strain rate field. The plots shown are for 
t ' = 0.4; in all cases the Argand number is 0. 

The degree to which the distribution of strain rates is controlled by the tem- 
perature dependence (Q /RT L ) is clearly seen. When Q /RT L is large (i.e. 50), 
high strain rates are confined to a region that moves progressively towards the 
interior of the fluid. For the intermediate value of Q /RTi (25), this separation 
is less pronounced or is absent. If it is present, the distance from the boundary of 
the locus of high strain rate is less than that for the equivalent experiment with 
larger Q /RT L . 

When Q I RT i is small, extensional strain rates are generally largest on the 
boundary of the box. In these instances, the change in Moho temperature 
required to affect the vertically averaged strength (see Eq. 4) is sufficiently large 
that it is not attained until near the end of the experiment (e.g. Pe — 300), or 
not attained at all (e.g., Pe = 1000 or 3000). 

The effect of lithospheric attenuation on the vertically averaged strength of 
the lithosphere and on the course of deformation is most clearly seen in the plots 
in which Q /RT L = 12.5 and the Peclet number is large; in this case a decrease 
in Moho temperature does not increase the strength. However, weakening of the 



- 15 - 


lithosphere due to thinning (see Eq. 3) causes strain to be concentrated in a small 
region along the boundary. Further extension there causes additional weakening; 
this positive feedback can be considered to be a kind of viscous stretching insta- 
bility, which for larger values of Q /RTi (> 25) is stabilized by the effect of 
temperature decreases brought about by thermal diffusion. 

Since the Peclet number is a measure of the relative rates of advection and 
diffusion of heat, it controls the timing of strength changes. Thus, it will affect 
when the locus of high strain rates begins to move away from the boundary, and 
how fast such a region will shift towards the interior of the fluid. This can also 
be seen in Figure 6. Experiments with small Peclet number, which have a fast 
rate of thermal diffusion relative to advection, show strain rate maxima that have 
had longer to develop and so are broader in area, shallower, and further from the 
boundary than experiments with higher Peclet numbers. In fact, for Pt = 300 
and Q /RT l = 25 or 50, the strain rate maximum has become so broad that it 
is not distinguishable with the contour interval chosen. 


More quantitatively, an understanding of the timing of strength changes can 
be obtained through a comparison of the results of the numerical calculation with 
the critical times predicted by the simple model of Section 3. The curves in Fig. 
7 are the critical times calculated as a fuction of strain rate with the right hand 
side of Eq. 7 evaluated for the parameters of Table 2, as discussed in Section 3. 
The relevant critical dimensionless time in the numerical experiments is chosen to 
be the time at which the strength B /Bq (Eq. 4) begins to increase at node 1 of 
the mesh (see Fig. 3 for node locations). The dimensionless strain rate in each 
case is taken to be the starting value at node 1 (« 2.2-see Fig. 4). These dimen- 
sionless times and strain rates are made dimensional for the purposes of plotting 


on Fig. 7 by multiplication or division by the appropriate timescale, which is 


D 2 
k Pe 


(see Secion 3). 
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The agreement between the critical times determined from the full thermo- 
mechanical calculations and those calculated from Eq. 10 is close enough to war- 
rant the use of Eq. 10 as a convenient rule-of-thumb. 

4.3.3. Vertical strain 

Figure 8 shows the calculated distribution of vertical strains (0) at a dimen- 
sionless time of 0.4 for experiments with different values of Pe and QfRT L . 
Substantial crustal thinning has occurred by this time for all experiments, but 
the distribution differs, owing to the effects of the temperature dependent rheol- 
ogy. 

Because of the migration of the locus of high strain rates away from the 
boundary (see Figs. 4 and 5) experiments with low Peclet numbers (in which 
rapid thermal diffusion has produced substantial horizontal differences in 
strength) tend to produce broad regions of moderately thinned crust. For exam- 
ple, when Pe — 300, strains are less than 0 = 2.5, except along the tapered part 
of the boundary. Conversely, when the Peclet number is large, the strain along 
the boundary exceeds 0 = 4 in all cases shown and the distribution is similar to 
that produced when the rheology is independent of temperature (Q /RT L = 0). 
To a lesser extent, these differences are also evident as a function of Q /RT L . A 
highly temperature dependent rheology results in more uniformly distributed 
extension than one with little or no temperature dependence. 

4.3.4. Buoyancy forces 

In the calculations presented so far, we have not considered the effects of 
buoyancy forces arising from horizontal variations in density. However, it is clear 
that such forces can be important whenever deformation produces lateral con- 
trasts in crustal or lithospheric thickness (e.g. Artyushkov, 1973; Fleitout and 
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Froidevaux, 1982; England and McKenzie, 1982; LePichon, 1983; Sonder et al., 
1987). During continental extension, thinning of the crust and attendant replace- 
ment by denser mantle gives rise to buoyancy forces that act to decrease the net 
potential energy available to drive continued extension. In contrast, if the man- 
tle lithosphere also thins, the replacement of mantle material with less dense 
asthenosphere has the opposite effect of tending to increase the extensional stress 
(see Le Pichon, 1983, Artyushkov, 1973). The net buoyancy force will be the sum 
of these two effects, and depends on the densities of crust, mantle, and astheno 
sphere, and the amount and vertical distribution of thinning. 

In terms of the thin viscous sheet model, the sensitivity of the lithosphere to 
such buoyancy forces is described by the Argand number (England and McKen- 
zie, 1982; also see the Appendix). Ar >0 corresponds to the situation in which 
buoyancy forces arising from horizontal variations in crustal thickness are large 
and opposite in sign to viscous forces driving the deformation. Such lithosphere 
can be considered too weak to support large gradients in crustal and lithospheric 
thickness. On the other hand, Ar <0 represents the state in which buoyancy 
forces resulting from extension act to augment other extensional forces. Finally, 
Ar =0 implies the state in which buoyancy forces do not influence the deforma- 
tion, neither damping nor amplifying the extension. 

Figure 9a shows the influence of the Argand number on the strain fields; 
experiments are shown for Ar — -3, 0, and 3, with Q /RTi = 0. Since when 
Ar > 0, buoyancy forces due to crustal thickness gradients result in pressure 
gradients directed away from areas of thicker crust and towards areas of thinner 
crust, less crustal thinning results than when Ar = 0. Conversely, when 
Ar < 0, the buoyancy forces act to localize deformation at previously thinned 
areas, so that deformation is confined to a smaller area in which strains are 
greater than when Ar >0. For example, when Ar = -3, the area with 0 > 4 



. - 18- 


includes approximately 2/3 of the total area with /? > 1.5, but only 1/3 when 
Ar = 0, and is confined to the vicinity of the taper in the boundary condition 
when Ar > 0. 

When Ar < 0 the deformation can be considered to be unstable, in the 
sense that it localizes around previously deformed regions. However, when a tem- 
perature dependent rheology is also included, the increase in strength of previ- 
ously strained areas stabilizes the deformation, as shown in Figure 9b, for experi- 
ments with Q JRT i = 50, and Ar =-3, 0, and 3. Although horizontal gra- 
dients in vertical strain are still greater when Ar < 0 (or less when Ar > 0) 
than when Ar = 0, in all cases the extensional strain and the horizontal gra- 
dients in strain at any point are less than when Q JRT i =0. 

5. Comparison with the Aegean 

We have presented general results for the extensional deformation of a thin 
viscous sheet with a temperature- and stress-dependent rheology. In this section 
we test the applicability of the model to continental extension by quantitatively 
comparing results of the calculations with observations of recent and active defor- 
mation in the Aegean region. 

Aegean extension began in the middle to late Miocene (13-10 Ma ago) (Le 
Pichon and Angelier, 1979) or possibly as late as 6 Ma ago (McKenzie, 1978b). 
Total extension of up to 100% is belived to have occurred in a north-south direc- 
tion (McKenzie, 1978b), with the greatest extensional strain located in the south- 
ern Aegean, north of Crete (e.g., Angelier et al., 1982). The region that has 
experienced extension (see Figure 11) is roughly 5X10 5 km 2 in area (Le Pichon, 
1983) and is bounded to the south by the Hellenic trench, and to the west by 
thrusting and compression in Albania and northwestern Greece. Extensional 
deformation dies out north of 42 'N in Yugoslavia, Bulgaria and northern 
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Greece, and east of 31 °E in Turkey (see Angelier et al., 1982). 

If the scale, D , is chosen to be 1000 km then the size of the extensional 
boundary in the calculations corresponds approximately to that of the Hellenic 
trench, and its displacement by 400 km at a dimensionless time of 0.4 is in rough 
agreement with the estimates quoted above for net extension in the Aegean. With 
this scaling, the range of Peclet numbers (300, 1000, 3000) investigated in the cal- 
culations of Section 4 correspond to magnitudes of 10, 30, and 100 mm yr -1 for 
the outward velocity of the boundary. The dimensionless time of 0.4, at which 
displacement of the extensional boundary reaches 400 km, corresponds to 40 Ma 
(Pe = 300), 13 Ma (Pe = 1000), and 4 Ma ( Pe — 3000). The longest of these 
times is greater than that believed to be available for the Aegean extension, but 
the two shorter times cover the range of estimates given by Le Pichon and 
Angelier (1979) and McKenzie (1978b) for the duration of the extension. 

The values of Q /RT L of 50, 25, and 12.5 correspond to activation energies 
for creep of 580, 290, and 150 kJ mol -1 . 

For convenience we shall refer to regions in the solution domain in terms of 
geographical directions — north, south, east, and west — with the +y -direction 
being north (see Figure 3). We emphasize that the geometry and boundary condi- 
tions shown in Figure 3 were chosen for their simplicity, and not to simulate in 
detail the kinematics of the Aegean. Hence we have not attempted to reproduce 
the exact geometry of the region, nor have we included any effects related to 
right lateral strike-slip motion along the North Anatolian fault. 

5.1. Extensional strain 

The pre-Miocene geology of the Aegean is similar to that of unextended 
regions in mainland Greece and Turkey (see summary by Horvath and Berckhe- 
mer, 1982). Therefore, rough estimates of the total extension in the Aegean can 
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be made by comparing present-day crustal thicknesses with those of Greece and 
Turkey, about 40-50 km (Makris, 1975; Makris and Vees, 1977). Seismic refrac- 
tion experiments have shown that the crustal thickness is 30-34 km under Crete, 
approximately 26 km beneath Amorgos in the central Aegean, and about 32 km 
under Evia, east of the Pelopenessos (Makris and Vees, 1977). This suggests thin- 
ning by factors of approximately 1.4, 1.7, and 1.4, respectively. 

Post-Miocene extensional strain in the Aegean has also been determined 
using plate motion data and subsidence data and is consistent with the rough 
estimates given above. From slip vectors of earthquakes located along the Hel- 
lenic trench, relative rotation rates of the European plate, African plate, and 
Aegean region and constraints imposed by fault trends and offsets, Le Pichon and 
Angelier (1979) estimated displacements along the boundaries of the extending 
region; the displacements were interpolated to give strains in the interior of the 
region. In the southern Aegean, average subsidence data for rectangular regions 
2100-3800 km 2 in area were used in conjunction with a simple uniform stretching 
model (McKenzie, 1978a) to obtain more detailed estimates of strain (Angelier et 
al., 1982). 

Uncertainties in the reconstructions are difficult to determine. Some indica- 
tion may be afforded by the calculated uncertainty in the strain averaged over 
the entire southern Aegean of /? = 1.65±0.25 given by Angelier et al (1982). 
However, it is more likely that strains are underestimated than overestimated. 
Determining the amount of extension from the ratio of the inferred pre- 
extensional crustal thickness to the present-day crustal thickness does not take 
into account additions to crustal thickness by mantle-derived magmatism, and 
determinations of strain that rely on offsets across normal faults will underesti- 
mate extension if the faults are listric. 
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To summarize the Aegean strain data, the region affected by extensional 
strain (/?>1.1) is approximately 700 km by 700 km in area. Maximum strain of 
f) = 1.8 occurs in the Sea of Crete; south of this a sharp decrease in strain (from 
f3 = 1.8 to 1.4) occurs over a distance about 100 km. Strains decrease gradually 
to the north of the Sea of Crete; the 0 = 1.1 contour is 400-500 km distant from 
the zone of maximum strain. Strains greater than f) = 1.5 are confined to the 
southern Aegean, and occur in an area approximately 250 km (N-S) by 500 km 
(E-W). 

To compare calculated strains (Fig. 8) with observed strains (Fig. 10), we 
consider calculated strains near the center of the extending region (i.e., along 
x — 0), since the large strains on the east and west margins of the area are 
numerical artefacts that occur due to the taper used in the boundary conditions. 
In all calculations, the portions of the fluid with 0 > 1.5 are within 450 km of 
the southern boundary, and in the case of the experiment with Pe = 300, 
Q /RT l = 50, and Ar — 0, the = 1.5 contour is about 250 km from the 
boundary. Maximum calculated strains along x' =0 range from less than 2.0 
at low Peclet number and high Q /RTl to considerably above 4.0 when the 
influence of temperature changes on the rheology is small (e.g. when Pe = 3000). 
None of the experimental results shows the sharp southwards decrease in from 
1.8 to 1.4 that occurs in the southern Aegean. The combination of parameters 
that exhibit a region with /? < 1.1 comparable in extent to that of the Aegean 
and maximum extensional strains /? < 2.5 is Pe ~ 1000, Q / RT L ^50 (c.f. Fig. 
8 ). 
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5.1.1. Strain rates 

A striking feature of the strain field in the Aegean is that while the greatest 
extensional strain is in the southernmost part of the region (Angelier et al., 1979; 
see Fig. 10, this paper) the most active region at present is some 300 km north of 
the trench. This is illustrated in Figure 11a, which shows the distribution of nor- 
mal and strike-slip earthquakes in the Aegean region for this century (Jackson 
and McKenzie, 1988), and Figure lib, which shows the shallow seismicity of all 
types for the region. 

Figure 12a shows the cumulative seismic moment release in the earthquakes 
of Fig. 11a, plotted as a function of distance north of 35 °N (the approximate 
southern limit of extensional deformation). Approximately 80% of the moment 
has been released in the region 400 km to 800 km north of the trench. As Figure 
12a and the analysis of Jackson and McKenzie (1988) shows, the average strain 
tensor expressed in this seismicity is nearly biaxial with north-south extension 
and vertical thinning; thus the profile of moment release in Fig. 12a may equally 
be regarded as a profile of (seismically expressed) northward displacement of the 
north of the region with respect to the south. Jackson and McKenzie (1988) esti- 
mate that the moment release shown in Fig. 12a corresponds to a displacement of 
2.4 to 8 m over the interval 1909-1981 or an average relative velocity of 33 to 113 
mm yr _1 . 

The observed profile of cumulative seismic moment release, or north-south 
displacement, may be compared with the velocity fields calculated in the experi- 
ments discussed above; this is done in Figures 13b, c, d. As in the case of the 
total strain, the quantity of interest is rather more sensitive to the Peclet number 
than to the value of Q /RT L . All calculations with Pe = 3000 show a concen- 
tration of strain rate towards the southern boundary, irrespective of the activa- 
tion energy shown. However, for each activation energy, a calculation with a 
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lower Peclet number can be found that exhibits slower strain rates in the south- 
ern 200-300 km of the extended region in the fashion shown by the seismic data; 
for Q /RTi = 50, this Peclet number is about 1000 (Fig. 12b), for 
Q jRT i = 25, it is between 1000 and 300 (Fig. 12c), and for Q /RT L = 12.5, it 
is about 300 (Fig. 12d). The best agreement between the normalized moment 
release curves and the calculated velocity curves is for Q /RT L ^ 50, 
Pe ~ 1000. 


5.2. Finite rotations and paleomagnetic data 

Figure 12 shows paleomagnetically determined rotations of Oligocene and 
younger rocks in the Aegean. Clockwise rotations of up to 25° are found for 
upper Miocene rocks on the Ionian islands of Kephallinia, Corfu, and Zakinthos 
(Laj et al., 1982). In the Epirus-Akarnania region of northwestern Greece, rota- 
tions of Eocene and Oligocene rocks range from 35-50 ° , clockwise (Horner and 
Freeman, 1983; Kissel et al., 1985). Clockwise rotations tend to be smaller east- 
wards: about 25° on the Chalkidiki peninsula (Kondopoulou and Westphal, 
1986), 26 ° on Skyros, 48 ° on Evia (Kissel et al., 1986a), and 6 ° on Lesbos 
(Kissel et al., 1986b). No significant rotation has occurred on Crete and Rhodes 
since the upper Miocene (Valente et al., 1982; Laj et al., 1982), in the Volos 
region in Greece since the middle Pliocene, or in Thrace since the lower Oligocene 
(Kissel et al., 1986a). Preliminary work on Miocene rocks in the Izmir area of 
western Turkey suggests 29 0 of counterclockwise rotation (Kissel et al., 1986b). 

Calculated finite rotations at t' = 0.4 are shown in Figure 13 for several 
numerical experiments. The results are not substantially different from each 
other, indicating that the finite rotations do not depend on the temperature 
dependence of the rheology, but are controlled mainly by the boundary condi- 
tions; the rotations of the material lines that make up the tapers on the 
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boundary account for most of the calculated rotations in the interior, with addi- 
tional rotation about vertical axes arising from the shear component of deforma- 
tion near these tapers. 

The calculated results agree with the paleomagnetic data: rotations are 
clockwise in the west, decrease to zero in the center of the extending region, and 
are counterclockwise in the east. However, the paleomagnetic data do not con- 
strain parameter ranges in the manner that the strain and strain rate fields do 
(Sections 5.1, 5.2). 

5.3. Summary of implications for the Aegean 

The combination of parameters that provide the closest agreement between 
the observed and calculated strains and strain rates is a Peclet number of around 
1000 and Q /RT L of 25 to 50. With the scaling of Table 2, the value of 1000 for 
the Peclet number corresponds to an outward velocity of 32 mm yr' 1 for the 
southern boundary, and the duration of the motion that produces displacement 
of this boundary by 400 km is 12.5 Ma. 

McKenzie (1987b) estimates 70 mm yr -1 for the velocity of southernmost 
Aegea with respect to Eurasia, and Jackson and McKenzie (1988) estimate values 
of 30-113 mm yr -1 from the seismic moment release in the region. The duration 
of 12.5 Ma for the extension is comparable with Le Pichon and Angelier’s (1979) 
estimate of 10-13 Ma and about twice McKenzie’s (1978b). Little significance 
should be attached, however, to a factor of two in these comparisons, owing to 
uncertainties in the observations and the simple configuration of the calculations. 

The average strain rate implied for the Aegean region from McKenzie’s 
(1978b) or Le Pichon and Angelier’s (1979) estimates of the amount and duration 
of strain is in the range 1 to 3 X1CT 15 s -1 . The scaling arguments of Section 3 
suggest that this is within the range of conditions in which thermal diffusion will 
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limit extension, and the rule-of-thumb illustrated in Fig. 2 shows that, for a litho- 
sphere whose rheology is determined by the creep of silicates, the maximum 
achievable extensional strain may be less than about 0 = 2. The comparisons of 
this section support the idea that the low strain rates and high strain in the 
southernmost Aegean are the result of this mechanism. 

6. Implications for stretching instabilities 

Recently, several workers (Fletcher and Hallet, 1983; Zuber and Parmentier, 
1986; Zuber et al., 1986) have studied lithospheric extension in terms of the 
mechanics of necking instabilities that develop in media consisting of horizontal 
layers with contrasting rheologies. In these models, the instabilities, once formed, 
grow exponentially with rates that depend on the strength contrasts of layers and 
with wavelengths that depend on the relative thicknesses of the layers. 

In this paper, we have identified two other phenomena that contribute to the 
instability of continental lithosphere during extension. First, the weakening of 
continental lithosphere in the early stages of extension (Fig. 4b), before there has 
been sufficient time for temperatures at or below the Moho to decrease, will tend 
to concentrate extension in previously thinned regions. Second, we have demon- 
strated that if buoyancy forces arising from horizontal density variations (Sec. 
4.3.4) act to increase the force available to drive extension ( Ar < 0), deforma- 
tion will be concentrated in narrower zones. 

However, the dominant effect of temperature dependence will be to stabilize 
continental extension. As the uppermost mantle temperatures decrease in 
stretched continental lithosphere, the strength will increase. It becomes more 
difficult to continue to extend already stretched lithosphere than to begin to 
extend unstretched lithosphere, and the growth of instabilities is highly damped. 
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7. Discussion and conclusions 

We have used a thin viscous sheet model to investigate the influence of a 
temperature-dependent rheology on continental deformation in an extensional 
enviroment. During the early stages of extension, attenuation of the lithosphere 
causes the vertically-averaged strength of the lithosphere to decrease slightly, and 
strain is concentrated into the regions that have already strained the most (e. g. 
Fig. 6; Fig. 8). However, after a critical time (which can be estimated from the 
arguments given in Section 3), thermal diffusion lowers the temperature in the 
lithosphere and so increases its strength. This critical time depends on the rela- 
tive rates of strain and thermal diffusion (described by the Peclet number) but is 
short compared with the thermal time constant of the lithosphere (Fig. 2). How- 
ever, at high strain rates the extensional strains that are achieved before this crit- 
ical time is reached are so great that it is reasonable to expect sea-floor spreading 
to begin. At lower strain rates, thermal diffusion can cause changes in strength at 
small extensional strains. These changes in strength cause the locus of maximum 
strain rate in a zone of extension to shift as much as a few hundred kilometers 
during the duration of extension; in general the location of the strain rate max- 
imum will move from regions of high strain to regions of lower strain. If 
sufficient time elapses, the effect of the migration of loci of high strain rate is the 
creation of wide zones of moderate extensional strain, rather than smaller areas of 
greater extension. The transition between essentially isothermal extension which, 
if continued, would lead to seafloor spreading, and extension that is almost com- 
pletely inhibited by strength increases due to thermal diffusion takes place 
between strain rates of 10 -14 and 10 -16 s -1 . 

Buoyancy forces due to horizontal gradients in crustal or lithospheric thick- 
ness can modify the deformation. Depending on whether the buoyancy forces act 
to augment (i.e. Ar < 0) or diminish (e.g. Ar > 0) the viscous forces arising 
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from the imposed boundary velocity, they can increase or decrease, respectively, 
the total strain and strain gradient at points within the deforming fluid. How- 
ever, when there is moderate or greater temperature dependence in the rheology 
(Q / RT i > 25), the effect of buoyancy forces on the deformation is of secondary 
importance. 

The Miocene-to recent extension of the Aegean region appears to have taken 
place at an average strain rate of 1 to 3X1CT 15 s -1 . Extensional strains are 
greatest in the southern Aegean, approaching 100% extension, but present-day 
strain rates are concentrated in the northern Aegean. This spatial separation of 
the regions of maximum strain and maximum strain rate may be the consequence 
of the temperature dependent rheology of lithospheric materials. The comparisons 
beween observations and calculations (Sections 5.1 and 5.2) reinforce this sugges- 
tion, and show that the present pattern of strain and strain rate in the Aegean 
are consistent with the deformation of lithosphere whose rheology is controlled by 
the properties of olivine determined by laboratory experiments (e.g. Goetze, 
1978), subjected to an extensional boundary moving at 30-60 mm yr~ l . 

In general, the effects of a temperature dependent rheology provide a 
mechanism for understanding why there can be spatial and temporal changes in 
the distribution of extension. This mechanism is a simple but inevitable conse- 
quence of the extension, and may reduce the need to invoke outside influences, 
such as changes in plate boundary forces, to explain observations of time- and 
space-varying continental extension. 
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Appendix: Physical and mathematical formulation 

This formulation follows that of England and McKenzie (1983) and House- 
man and England (1986a). The continental lithosphere is represented by a thin 
incompressible viscous fluid overlying an inviscid halfspace. Tractions are 
assumed to be zero on the top and bottom of the layer, which is in isostatic 
equilibrium. These assumptions allow the deformation of the sheet to be treated 
in terms of vertically averaged quantities. (Note that there is no consideration of 
the possibility of boudinage-type stretching instabilities of the lithosphere (see 
Zuber and Parmentier, 1986; Zuber et al., 1986). The deformation of the layer is 
governed by a vertically-averaged power law rheology (Eq. 2) where the tempera- 
ture dependence is included in the term B (Eq. 4). 

Given these assumptions, the equations of motion may be written in terms 
of horizontal components only: 


d 


dxi 



diij 

dxj 
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duj 

dii 



du dv 
dx dy 


^p c g(\-p c fp m ) 


dS 2 

dii 


(Al) 


where u and v are horizontal components of the velocity vector u, L is the 
lithosphere thickness, p c and p m are respectively the densities of crust and 
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mantle, S is the crustal thickness, g is the acceleration due to gravity, and i and 
j take on the values 1 and 2 only. The vertical component of velocity is obtained 
through the incompressibility condition (V‘U = 0). The effective viscosity, rj, is 


V— 




(A2) 


The variation of temperature, T , with position and time is described by the 
heat transfer equation 

4^- + U'VT = kv 2 T (A3) 

where k is thermal diffusivity, and t is time. Heat production is neglected. 
Equations (5.2) and (5.4) are made dimensionless by: 

(s',y',*',L',S' ) = [x,y,z,L,S)/L 0 
u' = u/v 0 
t 1 = t v q/L o 



T' = T/T l 

where B 0 is the value of B at t =0. The dimensionless equations are: 
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where primes have been omitted. In order to be consistent with England and 
McKenzie (1982), who used a slightly different nondimensionalization, we define 
the Argand number, Ar , as 

9 Pc (l-?t / 9m ) 

Ar = — (A8) 

BoK/Ao) " 


where L 0 is a characteristic vertical length scale, so that 



(AO) 


The Argand number may be thought of as the ratio of buoyancy forces arising 
from variations in crustal thickness of order L 0 to viscous forces required to 
deform the fluid at a reference strain rate v 0 /L 0 . When Ar = 0, forces due to 
crustal thickness contrasts do not affect the flow, and as Ar oo, the buoyancy 
forces dominate, making the fluid too weak to support large variations in crustal 
thickness, so that the deformation tends towards plane strain. As Ar — ► -oo, 
buoyancy forces also dominate, but they act in concert with the viscous forces so 
as to tend to enhance variations in crustal thickness. 


Pe , the Peclet number, equals , and is a ratio of the time constant for 

K 


thermal diffusion to the time constant for advection of heat. 


For a given crustal and lithospheric thickness distribution, Equation (A5) is 
solved, subject to the specified velocity or stress boundary conditions, to give the 
horizontal velocities, u and v . The vertical strain rate is then obtained through 
the incompressibility condition, and in a Lagrangean coordinate frame, the rate 
at which crustal thickness, S , increases is given by: 


1 dS • du dv 

S dt €zz dx dy , 


(A10) 
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A useful measure of strain is the finite rotation, which may be thought of as 
the rotation of an infinitetesimally small rigid disk embedded in the the deform- 
ing fluid (McKenzie and Jackson, 1983). It is obtained by integrating with respect 
to time half of the vertical component of the vorticity: 

»-/*«-£/( &-■£)« < A “) 


Numerical solution of governing equations 

The equations describing the mechanical behavior of the sheet (Eq. A5) are 
solved using a finite element method (see Houseman and England, 1986a, Appen- 
dix B). The solution domain is divided into triangular elements, as shown in Fig- 
ure 4. The Lagrangean coordinate system employed in the finite element 
approach permits considerable simplification of the calculation of the temperature 
field. In this reference frame the nodes travel with the fluid, so that there is no 
horizontal advection of heat relative to these nodes. In addition, horizontal gra- 
dients of temperature are small compared with vertical temperature gradients. 
Hence the dimensionless heat transfer equation reduces to: 

dT . _dT 1 d 2 T 

dt + * zzZ dz Pe dz 2 

At a given time t , the velocities in the thin layer are calculated from 
equation (A5) using the previously determined viscosities. The timestep At is 
then chosen so that the maximum strain in any element is less than or equal to 
2.5%. Convergence of the velocity solution has been checked by running calcula- 
tions with half the timestep or with approximately four times as many nodes. 
Once the timestep is chosen, a first order finite difference approximation to Eq. 
AIO allows the crustal thicknesses to be calculated at time t + At . 
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Temperatures are obtained by solving Eq. A6 using a Crank-Nicholson 
implicit finite difference scheme, which is uniformly stable. The finite difference 
grid spacing is L 0 /100. Convergence of the temperature calculation was checked 
by comparison with solutions obtained using half the number of grid points. 

After the new temperature field is obtained, the viscosity field at time t+/\t 
is found by numerically integrating Eq. 3 and substituting into Eq. A7. 
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Figure captions 

Figure 1: Sketch of thermal profiles of extended continental lithosphere, 
immediately after stretching (a) and a short time later (b). In the time interval 
between (a) and (b), thermal diffusion has lowered the Moho temperature by an 
amount AT, with a resultant increase in the strength of the lithosphere from B 0 
to B (see Eq. 4). The relevant length scale for the diffusion is the distance 
{ z L~ z m)IP between the base of the lithosphere and the base of the thinned crust. 
The original thicknesses of the lithosphere and crust are z L and z M . 

Figure 2: Time for cooling to influence strength of extending lithosphere 
plotted against strain rate, calculated from Eq. 10. Numbers on the curves refer 
to the constant in the denominator in the right hand side of Eq. 10. This con- 
stant depends on a number of parameters (see Eq. 7) but the range illustrated 
here probably covers most of the values appropriate to the extension of continen- 
tal lithosphere. 

Figure 3: Finite element mesh and boundary conditions. Velocities are cal- 
culated at the vertices and midpoints of the sides of the triangular elements. The 
nodes labelled 1--5 are those referred to in Section 4.3.1 and in Figure 4. The y - 
component of velocity on y = 0 is: v = -1.0 for j x | < 0.25 D ; 

v = -sin jn-(0.4-x)/0.3] for 0.252) < |x | <0.4 2); and v = 0 for 

0.4 2) < | x | . Initial crustal thickness is 0.035 X D. 

Figure 4: Moho temperature (a), strength (b), vertical strain rate (c), and 
vertical strain (d) for selected nodes, plotted as functions of time for an experi- 
ment with Pe = 300, Q /RTi = 50, Ar = 0, n =3. All units are dimension- 
less. Note that the vertical strain rate is plotted so that the absolute rate of 
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deformation increases upwards. Arrows on time axes correspond to the critical 
time for this run, calculated from Eq. 10. 

Figure 5: Plots of the distribution of dimensionless vertical strain rate ( e 22 ) 
at different times for the experiment shown in Figure 4 ( Pe — 300, 
Q / RT l == 50, Ar = 0, n = 3). 

Figure 6: Dimensionless vertical strain rate ( k 2i ) at t' = 0.4 for selected 
experiments. 

Figure 7: Comparison between the rule-of-thumb calculation of the time at 
which cooling of the lithosphere has increased its strength appreciably (Eq. 10) 
and the times at which extension rates begin to decrease in the calculations of 
Figs. 6 and 8. Parameter values are those of Table 2. Numbers on the curves 
refer to values of Q /RT L . 

Figure 8: Extensional strain (0) at t' = 0.4 for selected experiments. The 
Argand number is zero in all plots. Except for the /? = 1.1 contour, the contour 
interval is 0.5. 

Figure 9: Extensional strain (0) at t' = 0.4 for selected experiments with 
varying Argand number, (a) Q /RT L =0, Pe = 1000. (b) Q /RT L =50, 
Pe = 1000. 


Figure 10: Extensional strain estimated from observations in the Aegean 
(from Angelier et al., 1982). Contours indicate values of 0, the ratio of the thick- 
ness of undeformed crust to that of stretched crust. 
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Figure 11: Seismic activity in the Aegean, (a) Focal mechanisms of exten- 
sional and strike-slip earthquakes in the Aegean in the interval 1909-1981 (Jack- 
son and McKenzie, 1988). Larger mechanism represents the average of the indivi- 
dual mechanisms, each weighted by its scalar moment. 

(b) Epicenters of earthquakes in the Aegean whose hypocenters are shallower 
than 35 km and whose body wave magnitudes are greater than 4.5, from the 
USGS catalog over the interval 1966-1984. 

Figure 12: (a) Cumulative seismic moment for the events plotted in Fig. 
11a, as a function of distance north from 35 °N. Curve marked ”M 0 ” is for scalar 
moments and curve marked ”A/ yj is for the component of the moment tensor 
that represents extension in a north-south direction. The difference between the 
two curves primarily represents the contribution from strike-slip events within 
the region. 

(b-d) Comparison of the data of (a) with velocity profiles calculated along the y- 
axis of Fig. 3 for the calculations shown in Figs. 6 and 8. The seismic data from 
(a) are normalized so that the cumulative moment, whether Af 0 or M yy , is unity 
at the north of the region. Calculated velocity profiles are labelled with the Peclet 
number for the calculation and Q JRT L is 50 in (b), 25 in (c), and 12.5 in (d). 

Figure 13: Paleomagnetic rotations determined for locations in the Aegean. 
Wedges show the a 95 range about the mean direction. I, Ionian islands (Kefal- 
linia, Corfu, Zakinthos); E, Epirus-Akarnania; V, Volos; Ch,Chalkidiki peninsula; 
Ev, Evia; S, Skyros; T, Thrace; L, Lesbos; R, Rhodes; C, Crete; Iz, Izmir. 

Figure 14: Calculated finite rotations at an elapsed time of 13 Ma for 
selected experiments. Vectors show the amount of rotation, with the zero rotation 
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direction being north (towards the top of the page). Rotations less than 5° are 
not plotted. 



Table captions 


Table 1: Summary of dimensionless parameters and the ranges of values 
used in the calculations. 


Table 2: Values of constants used in the calculations. 



Table 1 


Parameter 

Description 

Value 

Pt 

Peclet number 

300, 1000, 3000 

Q/RT \ 

temperature dependence 

12.5, 25, 50 

At 

Argand number 

-3, 0, 3 

n 

stress exponent 

3 




Table 2 


Parameter Description 

Value 

D 

horizontal scale length 

1000 km 

L 0 

lithosphere thickness 

100 km 

So 

initial crustal thickness 

35 km 

*0 

time constant (= D Jv 0 ) 

33 Ma 

Pc 

density of crust 

3.0 X10 3 kg m -3 

Pm 

density of mantle 

3.3 X10 3 kg m -3 

9 

gravitational acceleration 

9.8 m s 2 

T l 

reference temperature 

1400 C 

K 

thermal diffusivity 

10 -6 m 2 s -1 

R 

gas constant 

g.SlXlO^kJ mol' 1 1C 1 






Depth 
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